This notebook is part of an exploratory analysis of machine learning used to decompose hyperspectral datasets of hybrid perovskite nanoscale materials. Two machine learning models are mainly used in this project: Nonnegative Matrix Factorization (NMF) and Variational Autoencoders (VAEs).
This notebook generates artificial hyperspectral data to test NMF and VAE models and evaluate their accuracy.
Notebook Six: Using Artificial Data
Imports, Functions, and Classes¶
Imports¶
from preprocessing import *
#import PIL
from PIL import ImageFilter
from PIL import Image
from scipy.special import wofz
import nimfa
import sklearn
from sklearn.decomposition import NMF
from sklearn import metrics
import torch; torch.manual_seed(0)
from torchvision.utils import make_grid
import torch.nn as nn
import torch.nn.functional as F
import torch.utils
import torch.distributions
import torchvision
import pyroved as pv
import pyro.distributions as dist
import pyro
from pyroved.models.base import baseVAE
from pyroved.nets import fcDecoderNet, fcEncoderNet, sDecoderNet
from pyroved.utils import (
generate_grid, get_sampler, set_deterministic_mode, to_onehot, transform_coordinates)
save_directory = "./svg-figures/artificial-data/" # modify to save figures as a specified file extension
file_extension = ".svg"
SEM images: f1_img1, f2_img1 CL images: f1_img2, f2_img2 Denoised data: f1_sb_median, f2_sb_median 2D denoised data: f1_denoised_2d, f2_denoised_2d Example points: f1_x_points, f1_y_points, f2_x_points, f2_y_points Wavelengths and dimensions: f1_wav, f2_wav, f1_xpix, f1_ypix, f1_zpix, f2_xpix, f2_ypix, f2_zpix
Functions¶
Distributions¶
# 3 distribution functions: gaussian, lorentzian, and voigt
def gaussian(x, amp1, cen1, sigma1):
return amp1 * (1 / (sigma1 * (np.sqrt(2 * np.pi)))) * (np.exp((-1.0 / 2.0) * (((x - cen1) / sigma1)**2)))
def lorentzian(x, amp1, cen1, wid1):
return (amp1 * wid1**2 / ((x - cen1)**2 + wid1**2))
def voigt(x, x0, y0, a, sig, gam):
return y0 + a * np.real(wofz((x - x0 + 1j * gam) / sig / np.sqrt(2))) / sig / np.sqrt(2 * np.pi)
NMF Functions¶
# SCIKIT-LEARN
def nmf_plot(suptitle, spect_matrices, img_matrices, xpix, ypix, save=False):
rows = len(spect_matrices)
columns = len(img_matrices)+3
gs = GridSpec(rows, columns)
gs.update(wspace=0.09,hspace=0)
figaspect = plt.figaspect(float(ypix*3)/float(xpix*6))
figsize = (5,8/(figaspect[0]/figaspect[1]))
fig = plt.figure(figsize=figsize)
#fig.patch.set_facecolor('#00000000')
fig.suptitle(suptitle, fontsize=10)
row_count = 0
for i in range(rows):
fig_plt = fig.add_subplot(gs[row_count,:2])
fig_plt.set_xlim([0, 299])
column_count = 0
for j in range(row_count + 2):
fig_img = fig.add_subplot(gs[row_count, column_count + 2])
fig_img.imshow(img_matrices[row_count][column_count, :].reshape(ypix, xpix))
fig_img.set_xticks([])
fig_img.set_yticks([])
fig_img.tick_params(color=color[j])
for spine in fig_img.spines.values():
spine.set_edgecolor(color[j])
spine.set(lw=2)
fig_plt.plot(spect_matrices[row_count][:, column_count], c=color[j], lw=1)
fig_plt.set_xticks([])
fig_plt.set_xticklabels([])
fig_plt.tick_params(axis='both', direction='out', length=3, width=1, labelsize=3)
if row_count == rows-1:
fig_plt.set_xticks([75, 150, 225])
fig_plt.set_xticklabels([75, 150, 225])
column_count += 1
row_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
plt.show()
def get_variance(model, data):
""" Estimate performance of the model on the data """
prediction = model.inverse_transform(model.transform(data))
return metrics.explained_variance_score(data, prediction, multioutput="variance_weighted")
def error_analysis(models, data):
temp = []
for model in models:
temp.append(get_variance(model.fit(data), data))
return temp
def explained_variance_plot(suptitle, rgb_results, hsi1_results, hsi_results, save=False):
fig = plt.figure(figsize=(4, 3))
#fig.patch.set_facecolor("#00000000")
gs = GridSpec(1, 1)
gs.update(wspace=0.1, hspace=0.1)
plt.rcParams["font.weight"] = "bold"
plt.rcParams["axes.labelweight"] = "bold"
fig_plt = fig.add_subplot(gs[0, 0])
fig_plt.plot([2, 3, 4, 5, 6], rgb_results[0], ms=10, marker='o', lw=2, c=color[0], label="Calibrator")
fig_plt.plot([2, 3, 4, 5, 6], hsi1_results[0], ms=10, marker='o', lw=2, c=color[1], label="Noiseless Data")
fig_plt.plot([2, 3, 4, 5, 6], hsi_results[0], ms=10, marker='o', lw=2, c=color[2], label="Noisy Data")
fig.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
fig_plt.set_title(suptitle)
fig_plt.set_ylabel("EVR")
fig_plt.set_xlabel("Number of Components")
fig_plt.set_ylim(0.98, 1.001)
fig_plt.set_xticks([2, 3, 4, 5, 6])
fig_plt.set_xticklabels([2, 3, 4, 5, 6])
fig_plt.tick_params(axis='both', direction='out', length=6, width=2, labelsize=6)
fig_plt.grid()
fig_plt.legend()
if save:
fig.savefig(save_directory + suptitle + file_extension, bbox_inches="tight")
# NIMFA
def nimfa_nmf_model(data, iterations, n_comp):
# Returns list of 2 matrices, W and H respectively
nmf = nimfa.Nmf(data, max_iter=iterations, rank=n_comp,
update=nimfa_nmf_update, objective='fro')
nmf_fit = nmf()
W = nmf_fit.basis()
H = nmf_fit.coef()
return [W, H]
def snmf_model(data, sparseness_level, n_comp, iterations, version):
# Returns list of 2 matrices, W and H respectively, and the snmf model
snmf = nimfa.Snmf(data, rank=n_comp, max_iter=iterations, version=version,
beta=sparseness_level)
snmf_fit = snmf()
W = snmf_fit.basis()
H = snmf_fit.coef()
return [W, H, snmf_fit.fit.sparseness()], snmf
def snmf_plot(suptitle, snmf_matrices, sparseness_levels, n_comp, xpix, ypix, save=False):
rows = len(sparseness_levels)
columns = n_comp + 2
gs = GridSpec(rows, columns)
gs.update(wspace=0.1,hspace=0.1)
figsize = (columns, rows)
fig = plt.figure(figsize=figsize)
fig.patch.set_facecolor('white')
fig.suptitle(suptitle, fontsize=10)
row_count = 0
for i in range(rows):
fig_plt = fig.add_subplot(gs[row_count,columns-2:])
fig_plt.set_xticks([])
fig_plt.set_yticks([])
fig_plt.tick_params(axis='both', direction='out', length=6, width=2)
column_count = 0
for j in range(columns-2):
fig_img = fig.add_subplot(gs[row_count,column_count])
fig_img.imshow(snmf_matrices[row_count][1][column_count,:].reshape(ypix, xpix))
fig_img.set_xticks([])
fig_img.set_yticks([])
fig_img.tick_params(color=color[j])
for spine in fig_img.spines.values():
spine.set_edgecolor(color[j])
spine.set(lw=2)
fig_plt.plot(snmf_matrices[row_count][0][:,column_count],
c=color[j], lw=1)
fig_plt.set_title(
"beta = {:} \nW sparseness = {:.3f} \nH sparseness = {:.3f}".format(
sparseness_levels[row_count],
snmf_matrices[row_count][2][0],
snmf_matrices[row_count][2][1]), y=1.0, pad=-20, fontsize=3)
column_count += 1
row_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
def error_analysis_nimfa(models):
temp = []
for model in models:
temp.append(model.evar())
return temp
def sparseness_evr_plot(suptitle, rgb_results, hsi1_results, hsi_results,
sparseness_levels, save=False):
points = len(sparseness_levels)
fig = plt.figure(figsize=(4, 3))
fig.suptitle(suptitle, fontsize=10)
#fig.patch.set_facecolor("#00000000")
gs = GridSpec(1, 1)
gs.update(wspace=0.1, hspace=0.1)
#plt.rcParams["font.weight"] = "bold"
#plt.rcParams["axes.labelweight"] = "bold"
fig_plt = fig.add_subplot(gs[0, 0])
fig_plt.plot(sparseness_levels, rgb_results[0], ms=10, marker='o', lw=2,
c=color[0], label="Calibrator")
fig_plt.plot(sparseness_levels, hsi1_results[0], ms=10, marker='o', lw=2,
c=color[1], label="Noiseless Data")
fig_plt.plot(sparseness_levels, hsi_results[0], ms=10, marker='o', lw=2,
c=color[2], label="Noisy Data")
fig.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
#fig_plt.set_title("Explained Variance for Blind NMF")
fig_plt.set_ylabel("EVR")
fig_plt.set_xlabel("Sparseness Level")
fig_plt.set_ylim(0.992, 1.001)
fig_plt.set_xscale("log")
fig_plt.set_xticks(sparseness_levels)
fig_plt.set_xticklabels(sparseness_levels)
fig_plt.tick_params(axis='both', direction='out', length=6, width=2, labelsize=6)
fig_plt.grid()
fig_plt.legend()
if save:
fig.savefig(save_directory + suptitle + file_extension, bbox_inches="tight")
plt.show()
VAE Functions¶
vae_plot builds a plot of the latent space dimensions with two or three latent dimensions and plots points corresponding to single pixels in the x and y axes of the original dataset. The point is to show a relationship between different spectral features detected by the autoencoder. Below the plot is an image component of the latent representation of the data.
spectra_plot plots each spectrum in the image color mapped according to its value along a specified axis.
ssvae_plot plots the latent space and image reconstruction of SSVAE models
pointfinder_plot shows every spectrum with a certain x value (along a vertical line in the image).
spectsig_plot shows all labeled spectra color mapped according to each class
ssvae_label_plot shows the first 180 spectra the SSVAE model assigns to a specified class
spect_label_plot plot one spectrum from each class in the labeled data
custom_spect_label_plot plot one spectrum from each class in the latent representation of the data
marker_ec = "black"
marker_lw = 0.08
marker_alpha = 0.5
marker_s = 10
# Plotting pixels containing entire spectra in a latent space and corresponding images
def vae_plot(suptitle, z_mean, xpix, ypix, save=False):
rows = 2
columns = len(z_mean[0])
fig = plt.figure(figsize=(4,4))
fig.suptitle(suptitle, fontsize=10)
fig.patch.set_facecolor("white")
for i in range(rows):
for j in range(columns):
if i == 0: # First row
if columns == 2: # 2D Latent space
ax = fig.add_subplot(2, 2, j + 1)
ax.scatter(z_mean[:, 0], z_mean[:, 1], c=z_mean[:, j], cmap=cmap,
ec=marker_ec, lw=marker_lw, alpha=marker_alpha, s=marker_s)
ax.tick_params(axis='both', which='major', labelsize=4, direction='out',
length=2, width=0.5)
else: # 3D Latent space
ax = fig.add_subplot(2, 3, j + 1, projection='3d')
ax.scatter(z_mean[:, 0], z_mean[:, 1], z_mean[:, 2],
c=z_mean[:, j], cmap=cmap,
ec=marker_ec, lw=marker_lw, alpha=marker_alpha, s=marker_s)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.zaxis.set_ticklabels([])
ax.set_xlabel("X", size=5, labelpad=-15)
ax.set_ylabel("Y", size=5, labelpad=-15)
ax.set_zlabel("Z", size=5, labelpad=-15)
if j == 0:
ax.view_init(30, -85)
elif j == 1:
ax.view_init(30, -40)
elif j == 2:
ax.view_init(30, 5)
if j == 0:
ax.set_title('Latent Space', size=7)
if i == 1: # Second row
if columns == 2: # 2D Latent Space
img = fig.add_subplot(2, 2, j+3)
else: # 3D Latent Space
img = fig.add_subplot(2, 3, j+4)
img.imshow(z_mean[:,j].reshape(ypix, xpix), cmap=cmap)
img.set_xticks([])
img.set_yticks([])
if save:
fig.savefig(save_directory + suptitle + file_extension)
# Spectra Plot (for determining latent dimensions)
def spectra_plot(suptitle, img, array, xpix, ypix, save=False):
n = xpix*ypix # total number of spectra
cmap = mpl.colormaps['viridis'](np.linspace(0,1,n)) # create colormap
fig = plt.figure(figsize=(6,3))
fig.suptitle(suptitle, fontsize=10)
fig.patch.set_facecolor('white')
ax = fig.add_subplot()
ax.set_prop_cycle('color', list(cmap)) # colormaps each spectra from first to last in array
ax.tick_params(axis='both', which='major', labelsize=5, direction='out', length=2, width=0.5)
#plt.xlim(short_wav - 5, long_wav + 5)
#ax.set_xticks([500, 750])
#ax.set_xticklabels([])
#plt.yscale('log')
for i in array:
ax.plot(img[:,i[0],i[1]], lw=0.03, alpha=0.5)
if save:
fig.savefig(save_directory + suptitle + file_extension)
# SEMI-SUPERVISED VAE
# Plot the latent space and image reconstruction of SSVAE models
def ssvae_plot(suptitle, z_mean, z_labels, xpix, ypix, save=False):
rows = 2
columns = 1
cmap = mpl.colormaps['viridis']
fig = plt.figure(figsize=(4,4))
fig.suptitle(suptitle, fontsize=10)
gs = GridSpec(rows, columns)
#fig.patch.set_facecolor('#00000000')
fig.patch.set_facecolor('white')
for i in range(rows):
column_count = 0
for j in range(columns):
if i == 0:
ax = fig.add_subplot(gs[0,column_count])
ax.scatter(z_mean[:,0], z_mean[:,1], c=z_labels, cmap=cmap,
ec=marker_ec, lw=marker_lw, alpha=marker_alpha, s=marker_s)
#ax.tick_params(axis='both', labelsize=10)
ax.set_xticks([])
ax.set_yticks([])
if i == 1:
img = fig.add_subplot(gs[1,column_count])
img.imshow(z_labels.reshape(ypix, xpix), cmap=cmap)
img.set_xticks([])
img.set_yticks([])
column_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
# Function to plot every spectrum with a certain x value (along a vertical line in the image)
def pointfinder_plot(img, x_pt, ypix):
rows = ypix
columns = 1
fig = plt.figure(figsize=(4, rows))
for i in range(rows):
ax = fig.add_subplot(rows, columns, i + 1)
ax.plot(img[:, i, x_pt], c=color[i], lw=2)
ax.text(0, 0, str([i, x_pt]), fontsize=4)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
# Plot all spectra in each of the 4 classes
def spectsig_plot(suptitle, img, pts, save=False):
cmap_list = mpl.colormaps['viridis'](np.linspace(0,1,len(pts))) # create colormap
rows = len(pts[0])
columns = 4
gs = GridSpec(rows, columns)
gs.update(wspace=0, hspace=0)
figsize = (columns*2, rows)
fig = plt.figure(figsize=figsize)
fig.patch.set_facecolor("white")
fig.suptitle(suptitle, fontsize=20)
#fig.subplots_adjust(top=0.95, bottom=0.1)
row_count = 0
for i in range(rows):
column_count = 0
for j in range(columns):
fig_plt = fig.add_subplot(gs[i,j])
fig_plt.plot(img[pts[j][i][0],pts[j]
[i][1],:], c=cmap_list[j], lw=2)
#plt.xlim(short_wav - 5, long_wav + 5) # Cropping the spectra shown in the plot
#fig_plt.tick_params(axis='both', direction='out', length=12, width=4)
fig_plt.set_xticks([])
fig_plt.set_yticks([])
column_count += 1
row_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
plt.show()
# This plot shows the first 180 spectra the SSVAE model assigns to a specified class
def ssvae_label_plot(img, pts, starting_index=0, save=False):
rows = 180
#rows = len(pts)
columns = 1
fig = plt.figure(figsize=(4, rows))
fig.patch.set_facecolor("white")
#fig.suptitle(suptitle, fontsize=30)
gs = GridSpec(rows, columns)
row_count = 0
for i in range(starting_index, starting_index + rows):
fig_plt = fig.add_subplot(gs[row_count,0])
fig_plt.plot(img[:, pts[i,0], pts[i,1]],
c=color[row_count], lw=3)
plt.xticks([])
plt.yticks([])
#plt.text(220,350,str((int(pts[i][0]),int(pts[i][1]))))
plt.text(0, 0, str(i))
row_count += 1
#gs.tight_layout(fig, rect=[0,0.03,1,0.95])
if save:
fig.savefig(save_directory + suptitle + file_extension)
plt.show()
# Plot one spectrum from each class in the labeled data
def spect_label_plot(suptitle, img, pts, save=False):
rows = 1
columns = 4
gs = GridSpec(rows, columns)
gs.update(wspace=0,hspace=0.2)
figsize = (columns*4, rows*3)
fig = plt.figure(figsize=figsize)
fig.suptitle(suptitle, size = 20)
fig.patch.set_facecolor("white")
#fig.patch.set_facecolor('#00000000')
#fig.suptitle(suptitle, fontsize=30)
#fig.subplots_adjust(top=0.95, bottom=0.1)
row_count = 0
for i in range(rows):
column_count = 0
for j in range(columns):
fig_plt = fig.add_subplot(gs[i,j])
fig_plt.plot(img[:,pts[j][i][0],pts[j][i][1]], c=cmap_list[j], lw=3)
#plt.text(300, 370, str(column_count))
#fig_plt.set_title(method_names[row_count] + point_names[column_count])
#fig_plt.tick_params(color=color[j])
#plt.xlim(short_wav - 5, long_wav + 5) # Cropping the spectra shown in the plot
#fig_plt.tick_params(axis='both', direction='out', length=12, width=4)
#fig_plt.set_xticks([400, 550, 750])
fig_plt.set_yticks([])
fig_plt.set_xticks([])
column_count += 1
row_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
plt.show()
# Plot one spectrum from each class in the latent representation of the data
def custom_spect_label_plot(suptitle, img, pts, save=False):
rows = 1
columns = 4
gs = GridSpec(rows, columns)
gs.update(wspace=0,hspace=0.2)
figsize = (columns*4, rows*3)
fig = plt.figure(figsize=figsize)
#fig.patch.set_facecolor("#00000000")
fig.patch.set_facecolor("white")
fig.suptitle(suptitle, fontsize=20)
#fig.subplots_adjust(top=0.95, bottom=0.1)
cmap_list = mpl.colormaps['viridis'](np.linspace(0,1,len(pts))) # create colormap
row_count = 0
for i in range(rows):
column_count = 0
for j in range(columns):
fig_plt = fig.add_subplot(gs[i, j])
fig_plt.plot(img[:, pts[j][0], pts[j][1]], c=cmap_list[j], lw=3)
#plt.text(300, 370, str(column_count))
#fig_plt.set_title(method_names[row_count] + point_names[column_count])
#fig_plt.tick_params(color=color[j])
#plt.xlim(short_wav, long_wav) # Cropping the spectra shown in the plot
fig_plt.tick_params(axis='both', direction='out', length=6, width=2)
#fig_plt.set_xticks([400, 750])
fig_plt.set_xticks([])
fig_plt.set_yticks([])
#fig_plt.set_xticklabels([])
#if column_count == 0:
#fig_plt.set_yticks([1000])
#fig_plt.set_yticklabels([])
#fig_plt.tick_params(axis='y', direction='out', length=12, width=4)
#if row_count == rows-1:
#fig_plt.set_xticks([400-short_wav, 750-long_wav])
#fig_plt.set_xticklabels([])
#fig_plt.tick_params(axis='x', direction='out', length=12, width=4)
column_count += 1
row_count += 1
if save:
fig.savefig(save_directory + suptitle + file_extension)
plt.show()
Hyperspectral Image Class¶
class HSI:
# Initialize the HSI
def __init__(self, img, dist1, dist2, dist3, zpix, ypix, xpix,
filter_type=None, filter_strength=0, noise=0):
self.filter_type = filter_type
self.filter_strength = filter_strength
# Add image filter
if self.filter_type == "Gaussian" or "gaussian":
self.im = img.filter(ImageFilter.GaussianBlur(radius = self.filter_strength))
if self.filter_type == "Box" or "box":
self.im = img.filter(ImageFilter.BoxBlur(radius = self.filter_strength))
if self.filter_type == None or "None" or "none":
self.im = img
# Initialize noise type
if noise == 0:
self.noise = 0
elif noise == 1:
self.noise = np.random.normal(0, 0.25, zpix)
elif noise == 2:
self.noise = np.random.normal(0, 0.5, zpix)
elif noise ==3 :
self.noise = np.random.normal(0, 1.0, zpix)
# Add noise to distributions
self.dist1 = dist1 + self.noise
self.dist2 = dist2 + self.noise
self.dist3 = dist3 + self.noise
# Create the HSI
r, g, b = self.im.split()
zero = r.point(lambda _ : 0) # create xpix by ypix array of zeros
red_merge = Image.merge("RGB", (r, zero, zero))
green_merge = Image.merge("RGB", (zero, g, zero))
blue_merge = Image.merge("RGB", (zero, zero, b))
red_array = np.asarray(red_merge)[:, :, 0]
green_array = np.asarray(green_merge)[:, :, 1]
blue_array = np.asarray(blue_merge)[:, :, 2]
red = normalize(red_array.flatten()).tolist()
blue = normalize(blue_array.flatten()).tolist()
green = normalize(green_array.flatten()).tolist()
for i in range(len(red)):
red[i] = red[i] * self.dist1
for i in range(len(green)):
green[i] = green[i] *self.dist2
for i in range(len(blue)):
blue[i] = blue[i] * self.dist3
hyper = []
for i in range(len(blue)):
idx = red[i] + green[i] + blue[i]
hyper.append(idx)
hyper = np.asarray(hyper).reshape(xpix, ypix, zpix)
hyper = np.swapaxes(hyper, 0, 2)
hyper = np.rot90(hyper, k=3, axes=(1, 2))
hyper = np.flip(hyper, axis=2)
self.hyper = hyper
# Original image
def img(self):
plt.imshow(self.im)
# Image at certain band
def hyper_img(self, band):
return(plt.imshow(self.hyper[:, :, band], cmap = "viridis"))
# Spectra at certain pixel
def hyper_band(self, x, y, color):
self.color = str(color)
fig, ax = plt.subplots(figsize=(2, 2))
ax.plot(self.hyper[:, y, x], color=self.color)
return fig
Generating Artificial Data¶
Reading Blazers Image¶
# Define image dimensions
zpix = 300
ypix = 150
xpix = 150
rgba = Image.open(r"png-figures/uab-blazers.png")
rgba = rgba.resize((xpix, ypix), Image.Resampling.LANCZOS)
#rgba = rgba.transpose(method=Image.FLIP_LEFT_RIGHT)
rgb = rgba.convert('RGB')
gray = rgb.convert('L')
# Display image
fig, ax = plt.subplots(figsize=(2, 2))
ax.imshow(rgb)
plt.tight_layout()
Create HSI¶
# distributions
x = np.linspace(0, zpix, zpix)
red_gaussian = gaussian(x, 1000, 75, 20)
green_lorentzian = lorentzian(x, 70, 150, 10)
blue_voigt = voigt(x, 225, 0, 500, 7, 1)
# HSI without noise
hsi1 = HSI(rgb, red_gaussian, green_lorentzian, blue_voigt, zpix, ypix, xpix)
hsi1 = hsi1.hyper
# HSI with Gaussian noise
hsi = HSI(rgb, red_gaussian, green_lorentzian, blue_voigt, zpix, ypix, xpix,
filter_type=None, filter_strength=0, noise=1)
hsi = hsi.hyper
#hsi.hyper_band(100, 200, color[0])
#plt.show()
fig, ax = plt.subplots(figsize=(2, 2))
ax.imshow(hsi[zpix//2, :, :])
plt.tight_layout()
Plot HSI¶
Plot Code¶
def artificial_data_plot(suptitle, hsi, pts, save=False):
# Create figure
fig, axes = plt.subplots(figsize=(10, 12))
fig.suptitle(suptitle, fontsize=20)
rows = 3
columns = 2
gs = GridSpec(rows, columns)
gs.update(wspace=0.1,hspace=0.1)
title_fontsize = 7
axes.spines["top"].set_visible(False)
axes.spines["right"].set_visible(False)
axes.spines["bottom"].set_visible(False)
axes.spines["left"].set_visible(False)
axes.set_xticks([])
axes.set_yticks([])
count = 0
for i in range(rows):
# img
ax = fig.add_subplot(gs[i, 0])
ax.imshow(hsi[pts[i][0][0], :, ], cmap=cmap)
ax.set_title("Wavelength: " + str(pts[i][0][0]), fontsize=title_fontsize)
ax.plot(pts[i][0][2], pts[i][0][1], "X", markersize=5, c=color[count])
ax.plot(pts[i][1][2], pts[i][1][1], "X", markersize=5, c=color[count + 1])
ax.set_xticks([])
ax.set_yticks([])
# plt
ax = fig.add_subplot(gs[i, 1])
ax.plot(hsi[:, pts[i][0][1], pts[i][0][2]], lw=2, c=color[count])
ax.plot(hsi[:, pts[i][1][1], pts[i][1][2]], lw=2, c=color[count + 1])
ax.axvline(x=pts[i][0][0], ls="--", c="black", lw=2)
count += 2
if save:
fig.savefig(save_directory + suptitle + file_extension)
Plots¶
# list of 3 sets of red and blue points ordered [z, y, x]
pts = [[[50, 75, 75], [50, 72, 110]], [[150, 22, 140], [150, 22, 60]],
[[250, 70, 20], [250, 70, 30]]]
artificial_data_plot("Artificial HSI with no Noise", hsi1, pts)
artificial_data_plot("Artificial HSI with Noise", hsi, pts)
Test¶
test = rgba.resize((150, 150), Image.Resampling.LANCZOS)
test = test.convert('RGB')
fig, ax = plt.subplots(figsize=(2, 2))
ax.imshow(test)
plt.tight_layout()
Extra¶
# Just in case
def artificial_data_plot(suptitle, hsi,
xr1=150, yr1=150, xb1=225, yb1=150, z1=230,
xr2=150, yr2=237, xb2=122, yb2=45, z2=420,
xr3=50, yr3=150, xb3=235, yb3=150, z3=800,
save=False):
# Create figure
fig, ax = plt.subplots(figsize=(10, 8))
fig.suptitle(suptitle, fontsize=20)
rows = 3
columns = 2
title_fontsize = 7
# img1
img1 = fig.add_subplot(rows, columns, 1)
img1.imshow(hsi.hyper[z1, :, ], cmap=cmap)
img1.set_title("Wavelength: " + str(z1), fontsize=title_fontsize)
img1.plot(xr1, yr1, "X", markersize=5, c=color[2]) # red
img1.plot(xb1, yb1, "X", markersize=5, c=color[0]) # blue
# plt1
plt1 = fig.add_subplot(rows, columns, 2)
plt1.plot(hsi.hyper[:, yr1, xr1], c=color[2]) # red
plt1.plot(hsi.hyper[:, yb1, xb1], c=color[0]) # blue
plt1.axvline(x=z1, ls="--", c=color[4], lw=2)
# img2
img2 = fig.add_subplot(rows, columns, 3)
img2.imshow(hsi.hyper[z2, :, ], cmap=cmap)
img2.set_title("Wavelength: " + str(z2), fontsize=title_fontsize)
img2.plot(xr2, yr2, "X", markersize=5, c=color[2]) # red
img2.plot(xb2, yb2, "X", markersize=5, c=color[0]) # blue
# plt2
plt2 = fig.add_subplot(rows, columns, 4)
plt2.plot(hsi.hyper[:, yr2, xr2], c=color[2]) # red
plt2.plot(hsi.hyper[:, yb2, xb2], c=color[0]) # blue
plt2.axvline(x=z2, ls="--", c=color[4], lw=2)
# img3
img3 = fig.add_subplot(rows, columns, 5)
img3.imshow(hsi.hyper[z3, :, ], cmap=cmap)
img3.set_title("Wavelength: " + str(z3), fontsize=title_fontsize)
img3.plot(xr3, yr3, "X", markersize=5, c=color[2]) # red
img3.plot(xb3, yb3, "X", markersize=5, c=color[0]) # blue
# plt3
plt3 = fig.add_subplot(rows, columns, 6)
plt3.plot(hsi.hyper[:, yr3, xr3], c=color[2]) # red
plt3.plot(hsi.hyper[:, yb3, xb3], c=color[0]) # blue
plt3.axvline(x=z3, ls="--", c=color[4], lw=2)
if save:
fig.savefig(save_directory + suptitle + file_extension)
Calibrator¶
# Read image
rgb_calibrator = Image.open(r"png-figures/rgb-calibrator.png")
rgb_calibrator = rgb_calibrator.resize((xpix, ypix),Image.Resampling.LANCZOS)
gray_calibrator = rgb_calibrator.convert('L')
# Display image
fig, ax = plt.subplots(figsize=(2, 2))
ax.imshow(rgb_calibrator)
plt.tight_layout()
# Create calibrator HSI
hsi_rgb = HSI(rgb_calibrator, red_gaussian, green_lorentzian, blue_voigt, zpix, ypix, xpix, noise=1)
hsi_rgb = hsi_rgb.hyper
fig, ax = plt.subplots(figsize=(2, 2))
ax.imshow(hsi_rgb[150, :, :])
plt.tight_layout()
calibrator_pts = [[[50, 75, 25], [50, 75, 75]], [[150, 75, 75], [150, 75, 125]],
[[250, 75, 25], [250, 75, 125]]]
artificial_data_plot("Calibrator HSI", hsi_rgb, calibrator_pts)
Preprocessing¶
Calibrator¶
hsi_rgb_1d = np.reshape(hsi_rgb, (zpix * ypix * xpix)) # reshape to 1D
# subtract background
# take mean of region with values close to or equal to zero
hsi_rgb_sb_mean = np.mean(hsi_rgb[0:49, 75, 125])
hsi_rgb_1d = [(i - hsi_rgb_sb_mean) for i in hsi_rgb_1d]
hsi_rgb_1d = [i.clip(min=0) for i in hsi_rgb_1d]
# add median filter
hsi_rgb_1d = median_filter(hsi_rgb_1d, size=3, footprint=None, output=None, mode="reflect", cval=0.0, origin=0)
hsi_rgb_denoised = np.reshape(hsi_rgb_1d, (zpix, ypix, xpix))
hsi_rgb_denoised_2d = np.reshape(hsi_rgb_1d, (zpix, ypix * xpix))
print(hsi_rgb_denoised[0:99, 75, 125])
[0. 0. 0. 0.27552875 0. 0.04978867 0. 0.16889373 0. 0. 0.38135244 0.20566114 0.10594008 0. 0. 0.09642526 0. 0. 0.10770786 0.19562126 0. 0.36177429 0. 0. 0.0503009 0.08127209 0.02828932 0. 0. 0.42239525 0.22861897 0.35820817 0.04569547 0. 0. 0.16142889 0.24382526 0. 0.28222269 0.0043404 0.07616022 0.15895464 0.01760716 0.13643076 0. 0. 0. 0.16622664 0.07762084 0.83714537 0. 0. 0. 0.04354635 0.43763247 0. 0.07563135 0. 0.09497708 0. 0. 0.32586182 0.33358066 0. 0. 0.75859235 0.17609968 0.06178614 0. 0.35788081 0.19915186 0.45181453 0.15375299 0. 0. 0.00478558 0. 0.1310737 0.01164888 0. 0.07062049 0.21185933 0.2387573 0. 0. 0. 0. 0.41862053 0.13288304 0.22579255 0.33791578 0. 0. 0. 0. 0. 0.04344398 0.08056229 0.29355606]
artificial_data_plot("Denoised Artificial HSI Calibrator", hsi_rgb_denoised, calibrator_pts)
Noisy data¶
hsi_1d = np.reshape(hsi, (zpix * ypix * xpix)) # reshape to 1D
# subtract background
hsi_sb_mean = np.mean(hsi[0:49, 70, 30])
hsi_1d = [(i - hsi_sb_mean) for i in hsi_1d]
hsi_1d = [i.clip(min=0) for i in hsi_1d]
# add median filter
hsi_1d = median_filter(hsi_1d, size=3, footprint=None, output=None, mode="reflect", cval=0.0, origin=0)
hsi_denoised = np.reshape(hsi_1d, (zpix, ypix, xpix))
hsi_denoised_2d = np.reshape(hsi_1d, (zpix, ypix * xpix))
print(hsi_denoised[0:49, 70, 30])
[0. 0. 0.08978075 0.00598767 0.09004533 0. 0. 0. 0. 0.16527938 0. 0.26038743 0.03756016 0. 0.11066113 0. 0. 0. 0.12138748 0.12681052 0. 0.03256904 0.12539103 0. 0. 0. 0. 0.06234503 0.28816074 0.29807671 0.13887105 0. 0.19668559 0.016202 0.29717581 0. 0. 0. 0.06047049 0.07211239 0.0576191 0.10714568 0. 0.19122916 0.1980068 0.04121401 0. 0. 0.0304507 ]
artificial_data_plot("Denoised Artificial HSI", hsi_denoised, pts)
NMF on Artificial Data¶
Blind NMF with Nimfa¶
Calibrator¶
iterations = 400
update = "divergence"
objective = "fro" # frobenius norm objective function
test_conv = 10 # how often convergence test is done
min_residuals = None # minimum required improvement of residuals from previous iteration
nmf0_2comp = nimfa.Nmf(hsi_rgb_denoised_2d, max_iter=iterations, rank=2, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf0_2comp_fit = nmf0_2comp()
nmf0_2comp_W = nmf0_2comp_fit.basis()
nmf0_2comp_H = nmf0_2comp_fit.coef()
nmf0_3comp = nimfa.Nmf(hsi_rgb_denoised_2d, max_iter=iterations, rank=3, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf0_3comp_fit = nmf0_3comp()
nmf0_3comp_W = nmf0_3comp_fit.basis()
nmf0_3comp_H = nmf0_3comp_fit.coef()
nmf0_4comp = nimfa.Nmf(hsi_rgb_denoised_2d, max_iter=iterations, rank=4, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf0_4comp_fit = nmf0_4comp()
nmf0_4comp_W = nmf0_4comp_fit.basis()
nmf0_4comp_H = nmf0_4comp_fit.coef()
nmf0_5comp = nimfa.Nmf(hsi_rgb_denoised_2d, max_iter=iterations, rank=5, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf0_5comp_fit = nmf0_5comp()
nmf0_5comp_W = nmf0_5comp_fit.basis()
nmf0_5comp_H = nmf0_5comp_fit.coef()
nmf0_6comp = nimfa.Nmf(hsi_rgb_denoised_2d, max_iter=iterations, rank=6, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf0_6comp_fit = nmf0_6comp()
nmf0_6comp_W = nmf0_6comp_fit.basis()
nmf0_6comp_H = nmf0_6comp_fit.coef()
# Arrays of matrices
nmf0_spect_matrices = [nmf0_2comp_W, nmf0_3comp_W, nmf0_4comp_W, nmf0_5comp_W, nmf0_6comp_W]
nmf0_img_matrices = [nmf0_2comp_H, nmf0_3comp_H, nmf0_4comp_H, nmf0_5comp_H, nmf0_6comp_H]
nmf_plot("Blind NMF on Calibrator 2-6 Components Nimfa", nmf0_spect_matrices,
nmf0_img_matrices, xpix, ypix)
Noiseless artificial data¶
hsi1_2d = np.reshape(hsi1, (zpix, ypix * xpix))
iterations = 400
update = "divergence"
objective = "fro" # frobenius norm objective function
test_conv = 10 # how often convergence test is done
min_residuals = None # minimum required improvement of residuals from previous iteration
nmf1_2comp = nimfa.Nmf(hsi1_2d, max_iter=iterations, rank=2, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf1_2comp_fit = nmf1_2comp()
nmf1_2comp_W = nmf1_2comp_fit.basis()
nmf1_2comp_H = nmf1_2comp_fit.coef()
nmf1_3comp = nimfa.Nmf(hsi1_2d, max_iter=iterations, rank=3, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf1_3comp_fit = nmf1_3comp()
nmf1_3comp_W = nmf1_3comp_fit.basis()
nmf1_3comp_H = nmf1_3comp_fit.coef()
nmf1_4comp = nimfa.Nmf(hsi1_2d, max_iter=iterations, rank=4, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf1_4comp_fit = nmf1_4comp()
nmf1_4comp_W = nmf1_4comp_fit.basis()
nmf1_4comp_H = nmf1_4comp_fit.coef()
nmf1_5comp = nimfa.Nmf(hsi1_2d, max_iter=iterations, rank=5, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf1_5comp_fit = nmf1_5comp()
nmf1_5comp_W = nmf1_5comp_fit.basis()
nmf1_5comp_H = nmf1_5comp_fit.coef()
nmf1_6comp = nimfa.Nmf(hsi1_2d, max_iter=iterations, rank=6, update=update, objective=objective,
test_conv=test_conv, min_residuals=min_residuals)
nmf1_6comp_fit = nmf1_6comp()
nmf1_6comp_W = nmf1_6comp_fit.basis()
nmf1_6comp_H = nmf1_6comp_fit.coef()
# Arrays of matrices
nmf1_spect_matrices = [nmf1_2comp_W, nmf1_3comp_W, nmf1_4comp_W, nmf1_5comp_W, nmf1_6comp_W]
nmf1_img_matrices = [nmf1_2comp_H, nmf1_3comp_H, nmf1_4comp_H, nmf1_5comp_H, nmf1_6comp_H]
#nmf1_spect_matrices = [nmf1_2comp_W]
#nmf1_img_matrices = [nmf1_2comp_H]
nmf_plot("Blind NMF on Noiseless Artificial Data 2-6 Components Nimfa", nmf1_spect_matrices,
nmf1_img_matrices, xpix, ypix)
Noisy artificial data¶
iterations = 400
update = "divergence"
objective = "fro" # frobenius norm objective function
test_conv = 10
min_residuals = None
nmf2_2comp = nimfa.Nmf(hsi_denoised_2d, max_iter=iterations, rank=2, update=update,
objective=objective, test_conv=test_conv, min_residuals=min_residuals)
nmf2_2comp_fit = nmf2_2comp()
nmf2_2comp_W = nmf2_2comp_fit.basis()
nmf2_2comp_H = nmf2_2comp_fit.coef()
nmf2_3comp = nimfa.Nmf(hsi_denoised_2d, max_iter=iterations, rank=3, update=update,
objective=objective, test_conv=test_conv, min_residuals=min_residuals)
nmf2_3comp_fit = nmf2_3comp()
nmf2_3comp_W = nmf2_3comp_fit.basis()
nmf2_3comp_H = nmf2_3comp_fit.coef()
nmf2_4comp = nimfa.Nmf(hsi_denoised_2d, max_iter=iterations, rank=4, update=update,
objective=objective, test_conv=test_conv, min_residuals=min_residuals)
nmf2_4comp_fit = nmf2_4comp()
nmf2_4comp_W = nmf2_4comp_fit.basis()
nmf2_4comp_H = nmf2_4comp_fit.coef()
nmf2_5comp = nimfa.Nmf(hsi_denoised_2d, max_iter=iterations, rank=5, update=update,
objective=objective, test_conv=test_conv, min_residuals=min_residuals)
nmf2_5comp_fit = nmf2_5comp()
nmf2_5comp_W = nmf2_5comp_fit.basis()
nmf2_5comp_H = nmf2_5comp_fit.coef()
nmf2_6comp = nimfa.Nmf(hsi_denoised_2d, max_iter=iterations, rank=6, update=update,
objective=objective, test_conv=test_conv, min_residuals=min_residuals)
nmf2_6comp_fit = nmf2_6comp()
nmf2_6comp_W = nmf2_6comp_fit.basis()
nmf2_6comp_H = nmf2_6comp_fit.coef()
# Arrays of matrices
nmf2_spect_matrices = [nmf2_2comp_W, nmf2_3comp_W, nmf2_4comp_W, nmf2_5comp_W, nmf2_6comp_W]
nmf2_img_matrices = [nmf2_2comp_H, nmf2_3comp_H, nmf2_4comp_H, nmf2_5comp_H, nmf2_6comp_H]
nmf_plot("Blind NMF on Noisy Artificial Data 2-6 Components Nimfa",
nmf2_spect_matrices, nmf2_img_matrices, xpix, ypix)
Blind NMF with Scikit-learn¶
init = "random"
solver = "cd"
beta_loss = "frobenius"
tolerance = 0
iterations = 400
random_state = None
alpha_W = 0.0
alpha_H = "same"
l1_ratio = 0.0
verbose = 0
shuffle = False
# 2 Components
nmf_2comp = NMF(n_components=2, init=init, solver=solver, beta_loss=beta_loss,
tol=tolerance, max_iter=iterations, random_state=random_state, alpha_W=alpha_W,
alpha_H=alpha_H, l1_ratio=l1_ratio, verbose=verbose, shuffle=shuffle)
# 3 Components
nmf_3comp = NMF(n_components=3, init=init, solver=solver, beta_loss=beta_loss,
tol=tolerance, max_iter=iterations, random_state=random_state, alpha_W=alpha_W,
alpha_H=alpha_H, l1_ratio=l1_ratio, verbose=verbose, shuffle=shuffle)
# 4 Components
nmf_4comp = NMF(n_components=4, init=init, solver=solver, beta_loss=beta_loss,
tol=tolerance, max_iter=iterations, random_state=random_state, alpha_W=alpha_W,
alpha_H=alpha_H, l1_ratio=l1_ratio, verbose=verbose, shuffle=shuffle)
# 5 Components
nmf_5comp = NMF(n_components=5, init=init, solver=solver, beta_loss=beta_loss,
tol=tolerance, max_iter=iterations, random_state=random_state, alpha_W=alpha_W,
alpha_H=alpha_H, l1_ratio=l1_ratio, verbose=verbose, shuffle=shuffle)
# 6 Components
nmf_6comp = NMF(n_components=6, init=init, solver=solver, beta_loss=beta_loss,
tol=tolerance, max_iter=iterations, random_state=random_state, alpha_W=alpha_W,
alpha_H=alpha_H, l1_ratio=l1_ratio, verbose=verbose, shuffle=shuffle)
Calibrator¶
# 2 Components
rgb_nmf_2comp = nmf_2comp
W_rgb_nmf_2comp = rgb_nmf_2comp.fit_transform(hsi_rgb_denoised_2d)
H_rgb_nmf_2comp = rgb_nmf_2comp.components_
# 3 Components
rgb_nmf_3comp = nmf_3comp
W_rgb_nmf_3comp = rgb_nmf_3comp.fit_transform(hsi_rgb_denoised_2d)
H_rgb_nmf_3comp = rgb_nmf_3comp.components_
# 4 Components
rgb_nmf_4comp = nmf_4comp
W_rgb_nmf_4comp = rgb_nmf_4comp.fit_transform(hsi_rgb_denoised_2d)
H_rgb_nmf_4comp = rgb_nmf_4comp.components_
# 5 Components
rgb_nmf_5comp = nmf_5comp
W_rgb_nmf_5comp = rgb_nmf_5comp.fit_transform(hsi_rgb_denoised_2d)
H_rgb_nmf_5comp = rgb_nmf_5comp.components_
# 6 Components
rgb_nmf_6comp = nmf_6comp
W_rgb_nmf_6comp = rgb_nmf_6comp.fit_transform(hsi_rgb_denoised_2d)
H_rgb_nmf_6comp = rgb_nmf_6comp.components_
# Spectral components
rgb_nmf_spect_matrices = [W_rgb_nmf_2comp, W_rgb_nmf_3comp, W_rgb_nmf_4comp, W_rgb_nmf_5comp,
W_rgb_nmf_6comp]
# Image components
rgb_nmf_img_matrices = [H_rgb_nmf_2comp, H_rgb_nmf_3comp, H_rgb_nmf_4comp, H_rgb_nmf_5comp,
H_rgb_nmf_6comp]
nmf_plot("Blind NMF on Calibrator 2-6 Components",
rgb_nmf_spect_matrices, rgb_nmf_img_matrices, xpix, ypix)
Noiseless artificial data¶
# 2 Components
hsi1_nmf_2comp = nmf_2comp
W_hsi1_nmf_2comp = hsi1_nmf_2comp.fit_transform(hsi1_2d)
H_hsi1_nmf_2comp = hsi1_nmf_2comp.components_
# 3 Components
hsi1_nmf_3comp = nmf_3comp
W_hsi1_nmf_3comp = hsi1_nmf_3comp.fit_transform(hsi1_2d)
H_hsi1_nmf_3comp = hsi1_nmf_3comp.components_
# 4 Components
hsi1_nmf_4comp = nmf_4comp
W_hsi1_nmf_4comp = hsi1_nmf_4comp.fit_transform(hsi1_2d)
H_hsi1_nmf_4comp = hsi1_nmf_4comp.components_
# 5 Components
hsi1_nmf_5comp = nmf_5comp
W_hsi1_nmf_5comp = hsi1_nmf_5comp.fit_transform(hsi1_2d)
H_hsi1_nmf_5comp = hsi1_nmf_5comp.components_
# 6 Components
hsi1_nmf_6comp = nmf_6comp
W_hsi1_nmf_6comp = hsi1_nmf_6comp.fit_transform(hsi1_2d)
H_hsi1_nmf_6comp = hsi1_nmf_6comp.components_
# Spectral components
hsi1_nmf_spect_matrices = [W_hsi1_nmf_2comp, W_hsi1_nmf_3comp, W_hsi1_nmf_4comp, W_hsi1_nmf_5comp,
W_hsi1_nmf_6comp]
# Image components
hsi1_nmf_img_matrices = [H_hsi1_nmf_2comp, H_hsi1_nmf_3comp, H_hsi1_nmf_4comp, H_hsi1_nmf_5comp,
H_hsi1_nmf_6comp]
nmf_plot("Blind NMF on Noiseless Artificial Data 2-6 Components",
hsi1_nmf_spect_matrices, hsi1_nmf_img_matrices, xpix, ypix)
Noisy artificial data¶
# 2 Components
hsi_nmf_2comp = nmf_2comp
W_hsi_nmf_2comp = hsi_nmf_2comp.fit_transform(hsi_denoised_2d)
H_hsi_nmf_2comp = hsi_nmf_2comp.components_
# 3 Components
hsi_nmf_3comp = nmf_3comp
W_hsi_nmf_3comp = hsi_nmf_3comp.fit_transform(hsi_denoised_2d)
H_hsi_nmf_3comp = hsi_nmf_3comp.components_
# 4 Components
hsi_nmf_4comp = nmf_4comp
W_hsi_nmf_4comp = hsi_nmf_4comp.fit_transform(hsi_denoised_2d)
H_hsi_nmf_4comp = hsi_nmf_4comp.components_
# 5 Components
hsi_nmf_5comp = nmf_5comp
W_hsi_nmf_5comp = hsi_nmf_5comp.fit_transform(hsi_denoised_2d)
H_hsi_nmf_5comp = hsi_nmf_5comp.components_
# 6 Components
hsi_nmf_6comp = nmf_6comp
W_hsi_nmf_6comp = hsi_nmf_6comp.fit_transform(hsi_denoised_2d)
H_hsi_nmf_6comp = hsi_nmf_6comp.components_
# Spectral components
hsi_nmf_spect_matrices = [W_hsi_nmf_2comp, W_hsi_nmf_3comp, W_hsi_nmf_4comp, W_hsi_nmf_5comp,
W_hsi_nmf_6comp]
# Image components
hsi_nmf_img_matrices = [H_hsi_nmf_2comp, H_hsi_nmf_3comp, H_hsi_nmf_4comp, H_hsi_nmf_5comp,
H_hsi_nmf_6comp]
nmf_plot("Blind NMF on Noisy Artificial Data 2-6 Components",
hsi_nmf_spect_matrices, hsi_nmf_img_matrices, xpix, ypix)
Explained Variance¶
rgb_nmf_models = [rgb_nmf_2comp, rgb_nmf_3comp, rgb_nmf_4comp, rgb_nmf_5comp, rgb_nmf_6comp]
rgb_nmf_results = [error_analysis(rgb_nmf_models, hsi_rgb_denoised_2d)]
hsi1_nmf_models = [hsi1_nmf_2comp, hsi1_nmf_3comp, hsi1_nmf_4comp, hsi1_nmf_5comp, hsi1_nmf_6comp]
hsi1_nmf_results = [error_analysis(hsi1_nmf_models, hsi1_2d)]
hsi_nmf_models = [hsi_nmf_2comp, hsi_nmf_3comp, hsi_nmf_4comp, hsi_nmf_5comp, hsi_nmf_6comp]
hsi_nmf_results = [error_analysis(hsi_nmf_models, hsi_denoised_2d)]
explained_variance_plot("Explained Variance of Blind NMF on Artificial Data", rgb_nmf_results,
hsi1_nmf_results, hsi_nmf_results)
Sparse NMF with Nimfa¶
n_comp = 3
iterations = 400
sparseness_levels = [1e-8, 1e-6, 1e-4, 1e-3, 0.01, 0.1, 1.0, 5.0]
rgb_snmf_matrices = []
hsi1_snmf_matrices = []
hsi_snmf_matrices = []
rgb_snmf_models = []
hsi1_snmf_models = []
hsi_snmf_models = []
for s in sparseness_levels: # all models enforce sparseness in the spectral components
rgb_matrices, rgb_model = (snmf_model(hsi_rgb_denoised_2d, s, n_comp, iterations, "l"))
rgb_snmf_matrices.append(rgb_matrices)
rgb_snmf_models.append(rgb_model)
print("Calibrator SNMF with sparseness level " + str(s) + " completed")
hsi1_matrices, hsi1_model = (snmf_model(hsi1_2d, s, n_comp, iterations, "l"))
hsi1_snmf_matrices.append(hsi1_matrices)
hsi1_snmf_models.append(hsi1_model)
print("Noiseless SNMF with sparseness level " + str(s) + " completed")
hsi_matrices, hsi_model = (snmf_model(hsi_denoised_2d, s, n_comp, iterations, "l"))
hsi_snmf_matrices.append(hsi_matrices)
hsi_snmf_models.append(hsi_model)
print("Noisy SNMF with sparseness level " + str(s) + " completed")
Calibrator SNMF with sparseness level 1e-08 completed Noiseless SNMF with sparseness level 1e-08 completed Noisy SNMF with sparseness level 1e-08 completed Calibrator SNMF with sparseness level 1e-06 completed Noiseless SNMF with sparseness level 1e-06 completed Noisy SNMF with sparseness level 1e-06 completed Calibrator SNMF with sparseness level 0.0001 completed Noiseless SNMF with sparseness level 0.0001 completed Noisy SNMF with sparseness level 0.0001 completed Calibrator SNMF with sparseness level 0.001 completed Noiseless SNMF with sparseness level 0.001 completed Noisy SNMF with sparseness level 0.001 completed Calibrator SNMF with sparseness level 0.01 completed Noiseless SNMF with sparseness level 0.01 completed Noisy SNMF with sparseness level 0.01 completed Calibrator SNMF with sparseness level 0.1 completed Noiseless SNMF with sparseness level 0.1 completed Noisy SNMF with sparseness level 0.1 completed Calibrator SNMF with sparseness level 1.0 completed Noiseless SNMF with sparseness level 1.0 completed Noisy SNMF with sparseness level 1.0 completed Calibrator SNMF with sparseness level 5.0 completed Noiseless SNMF with sparseness level 5.0 completed Noisy SNMF with sparseness level 5.0 completed
Calibrator¶
snmf_plot("Sparse NMF with Calibrator",
rgb_snmf_matrices, sparseness_levels, n_comp, xpix, ypix)
Noiseless artificial data¶
snmf_plot("Sparse NMF with Noiseless Artificial Data",
hsi1_snmf_matrices, sparseness_levels, n_comp, xpix, ypix)
Noisy artificial data¶
snmf_plot("Sparse NMF with Noisy Artificial Data",
hsi_snmf_matrices, sparseness_levels, n_comp, xpix, ypix)
Sparse NMF Explained Variance¶
rgb_snmf_results = [error_analysis_nimfa(rgb_snmf_models)]
hsi1_snmf_results = [error_analysis_nimfa(hsi1_snmf_models)]
hsi_snmf_results = [error_analysis_nimfa(hsi_snmf_models)]
sparseness_evr_plot("Explained Variance of Sparse NMF on Artificial Data", rgb_snmf_results,
hsi1_snmf_results, hsi_snmf_results, sparseness_levels)
Unsupervised VAE on Artificial Data¶
Calibrator¶
# Transpose and reshape
rgb_denoised_1d_transpose = hsi_rgb_denoised_2d.transpose(1,0).reshape(zpix * ypix * xpix)
rgb_denoised_2d_norm = normalize(rgb_denoised_1d_transpose).reshape(xpix * ypix, zpix)
# Convert to Torch tensor object
rgb_denoised_2d_norm_t = torch.from_numpy(np.array(rgb_denoised_2d_norm).astype('float64')).float()
rgb_n_samples = rgb_denoised_2d_norm_t.size()[0] # number of spectral points
rgb_l_signal = rgb_denoised_2d_norm_t.size()[1] # number of spectra
rgb_train_data = rgb_denoised_2d_norm_t.clone()
rgb_train_loader = pv.utils.init_dataloader(rgb_train_data.unsqueeze(1), batch_size=64)
rgb_in_dim = (rgb_l_signal,)
# Gradually increase the weight of the KL Divergence term in the loss function with each epoch
rgb_kl_scale = torch.linspace(start=.001, end=.01, steps=50)
# Train the model
rgb_vae = pv.models.iVAE(rgb_in_dim, latent_dim=2, invariances=None, sampler_d="gaussian")
rgb_vae_trainer = pv.trainers.SVItrainer(rgb_vae)
for e in range(50):
sc = rgb_kl_scale[e] if e < len(rgb_kl_scale) else rgb_kl_scale[-1]
rgb_vae_trainer.step(rgb_train_loader, scale_factor=sc)
if e//10 == e/10:
rgb_vae_trainer.print_statistics()
rgb_vae_z_mean, rgb_vae_z_sd = rgb_vae.encode(rgb_train_data)
rgb_vae_z_mean = normalize(rgb_vae_z_mean)
rgb_vae_z_sd = normalize(rgb_vae_z_sd)
Epoch: 1 Training loss: 74.4753 Epoch: 11 Training loss: 67.7587 Epoch: 21 Training loss: 67.7604 Epoch: 31 Training loss: 67.7679 Epoch: 41 Training loss: 67.7786
vae_plot("Vanilla VAE on Calibrator", rgb_vae_z_mean, xpix, ypix)
Noiseless artificial data¶
# Transpose and reshape
hsi1_denoised_1d_transpose = hsi1_2d.transpose(1,0).reshape(zpix * ypix * xpix)
hsi1_denoised_2d_norm = normalize(hsi1_denoised_1d_transpose).reshape(xpix * ypix, zpix)
# Convert to Torch tensor object
hsi1_denoised_2d_norm_t = torch.from_numpy(np.array(hsi1_denoised_2d_norm).astype('float64')).float()
hsi1_n_samples = hsi1_denoised_2d_norm_t.size()[0] # number of spectral points
hsi1_l_signal = hsi1_denoised_2d_norm_t.size()[1] # number of spectra
hsi1_train_data = hsi1_denoised_2d_norm_t.clone()
hsi1_train_loader = pv.utils.init_dataloader(hsi1_train_data.unsqueeze(1), batch_size=64)
hsi1_in_dim = (hsi1_l_signal,)
# Gradually increase the weight of the KL Divergence term in the loss function with each epoch
hsi1_kl_scale = torch.linspace(start=.001, end=.01, steps=50)
# Train the model
hsi1_vae = pv.models.iVAE(hsi1_in_dim, latent_dim=2, invariances=None, sampler_d="gaussian")
hsi1_vae_trainer = pv.trainers.SVItrainer(hsi1_vae)
for e in range(50):
sc = hsi1_kl_scale[e] if e < len(hsi1_kl_scale) else hsi1_kl_scale[-1]
hsi1_vae_trainer.step(hsi1_train_loader, scale_factor=sc)
if e//10 == e/10:
hsi1_vae_trainer.print_statistics()
hsi1_vae_z_mean, hsi1_vae_z_sd = hsi1_vae.encode(hsi1_train_data)
hsi1_vae_z_mean = normalize(hsi1_vae_z_mean)
hsi1_vae_z_sd = normalize(hsi1_vae_z_sd)
Epoch: 1 Training loss: 73.6469 Epoch: 11 Training loss: 67.7582 Epoch: 21 Training loss: 67.7607 Epoch: 31 Training loss: 67.7626 Epoch: 41 Training loss: 67.7671
vae_plot("Vanilla VAE on Noiseless Artificial Data", hsi1_vae_z_mean, xpix, ypix)
hsi1_vae.manifold2d(d=5)
tensor([[2.0481e-03, 2.0453e-03, 2.1140e-03, ..., 2.2354e-03, 2.1738e-03,
2.1642e-03],
[1.4051e-03, 1.4795e-03, 1.3848e-03, ..., 1.4651e-03, 1.5669e-03,
1.4010e-03],
[1.2660e-04, 1.2518e-04, 1.3661e-04, ..., 1.1144e-04, 1.1084e-04,
1.0773e-04],
...,
[1.0392e-04, 8.0720e-05, 9.9335e-05, ..., 9.0773e-05, 9.5632e-05,
8.1836e-05],
[5.1150e-05, 4.4734e-05, 5.3905e-05, ..., 5.1767e-05, 5.4062e-05,
5.0325e-05],
[4.9471e-05, 4.7185e-05, 5.2956e-05, ..., 5.4603e-05, 5.4376e-05,
5.2619e-05]])
Noisy artificial data¶
# Transpose and reshape
hsi_denoised_1d_transpose = hsi_denoised_2d.transpose(1,0).reshape(zpix * ypix * xpix)
hsi_denoised_2d_norm = normalize(hsi_denoised_1d_transpose).reshape(xpix * ypix, zpix)
# Convert to Torch tensor object
hsi_denoised_2d_norm_t = torch.from_numpy(np.array(hsi_denoised_2d_norm).astype('float64')).float()
hsi_n_samples = hsi_denoised_2d_norm_t.size()[0] # number of spectral points
hsi_l_signal = hsi_denoised_2d_norm_t.size()[1] # number of spectra
hsi_train_data = hsi_denoised_2d_norm_t.clone()
hsi_train_loader = pv.utils.init_dataloader(hsi_train_data.unsqueeze(1), batch_size=64)
hsi_in_dim = (hsi_l_signal,)
# Gradually increase the weight of the KL Divergence term in the loss function with each epoch
hsi_kl_scale = torch.linspace(start=.001, end=.01, steps=50)
# Train the model
hsi_vae = pv.models.iVAE(hsi_in_dim, latent_dim=2, invariances=None, sampler_d="gaussian")
hsi_vae_trainer = pv.trainers.SVItrainer(hsi_vae)
for e in range(50):
sc = hsi_kl_scale[e] if e < len(hsi_kl_scale) else hsi_kl_scale[-1]
hsi_vae_trainer.step(hsi_train_loader, scale_factor=sc)
if e//10 == e/10:
hsi_vae_trainer.print_statistics()
hsi_vae_z_mean, hsi_vae_z_sd = hsi_vae.encode(hsi_train_data)
hsi_vae_z_mean = normalize(hsi_vae_z_mean)
hsi_vae_z_sd = normalize(hsi_vae_z_sd)
Epoch: 1 Training loss: 73.6643 Epoch: 11 Training loss: 67.7575 Epoch: 21 Training loss: 67.7588 Epoch: 31 Training loss: 67.7622 Epoch: 41 Training loss: 67.7672
vae_plot("Vanilla VAE on Noisy Artificial Data", hsi_vae_z_mean, xpix, ypix)
Determining latent dimensions¶
# X-axis
# Convert to np array and reshape to 2d
hsi_vae_z_mean_x = hsi_vae_z_mean[:, 0].detach().cpu().numpy()
hsi_vae_z_mean_x = hsi_vae_z_mean_x.reshape(ypix, xpix)
print("X Min: ", (hsi_vae_z_mean_x==np.min(hsi_vae_z_mean_x)).nonzero())
print("X Max: ", (hsi_vae_z_mean_x==np.max(hsi_vae_z_mean_x)).nonzero())
hsi_vae_z_mean_x_sort = np.dstack(
np.unravel_index(np.argsort(hsi_vae_z_mean_x.ravel()), (ypix, xpix)))
hsi_vae_z_mean_x_sort = np.squeeze(hsi_vae_z_mean_x_sort)
spectra_plot("Vanilla VAE on Noisy Artificial Data X-axis Spectra", hsi, hsi_vae_z_mean_x_sort,
xpix, ypix, True)
X Min: (array([86]), array([112])) X Max: (array([85, 85]), array([81, 82]))
# Y-axis
hsi_vae_z_mean_y = hsi_vae_z_mean[:, 1].detach().cpu().numpy()
hsi_vae_z_mean_y = hsi_vae_z_mean_y.reshape(ypix, xpix)
print("Y Min: ", (hsi_vae_z_mean_y==np.min(hsi_vae_z_mean_x)).nonzero())
print("Y Max: ", (hsi_vae_z_mean_y==np.max(hsi_vae_z_mean_x)).nonzero())
hsi_vae_z_mean_y_sort = np.dstack(
np.unravel_index(np.argsort(hsi_vae_z_mean_y.ravel()), (ypix, xpix)))
hsi_vae_z_mean_y_sort = np.squeeze(hsi_vae_z_mean_y_sort)
spectra_plot("Vanilla VAE on Noisy Artificial Data Y-axis Spectra", hsi, hsi_vae_z_mean_y_sort,
xpix, ypix, True)
Y Min: (array([], dtype=int64), array([], dtype=int64)) Y Max: (array([], dtype=int64), array([], dtype=int64))
hsi_vae.manifold2d(d=5)
tensor([[7.1727e-05, 4.1761e-05, 1.2139e-03, ..., 2.9731e-03, 1.2540e-04,
3.7594e-05],
[6.7384e-05, 4.2683e-05, 8.1973e-04, ..., 1.9013e-03, 9.5944e-05,
3.8472e-05],
[6.2597e-05, 6.1596e-05, 1.5167e-04, ..., 1.3518e-04, 4.9536e-05,
5.5334e-05],
...,
[4.1632e-05, 4.9906e-05, 9.3308e-05, ..., 8.5248e-05, 3.2839e-05,
5.0167e-05],
[3.4235e-05, 3.5931e-05, 4.9513e-05, ..., 5.3711e-05, 2.1183e-05,
3.4411e-05],
[3.6095e-05, 3.1962e-05, 4.2673e-05, ..., 5.2135e-05, 1.7393e-05,
3.5447e-05]])
Semi-Supervised VAE on Artificial Data¶
Preparing data¶
pointfinder_plot(hsi, 55, ypix)
# arrays of points labeled 0 to 3, totalling 4 classes
# format: [y, x]
pts0 = [[0, 0], [0, 50], [0, 100], [94, 70], [115, 30], [127, 30]]
pts1 = [[6, 50], [8, 50], [51, 70], [69, 70], [50, 30], [58, 30]]
pts2 = [[17, 50], [19, 50], [104, 70], [113, 70], [62, 30], [70, 30]]
pts3 = [[56, 60], [47, 70], [49, 60], [53, 60], [47, 55], [53, 55]]
labels = [pts0, pts1, pts2, pts3]
hsi_denoised_3d_norm_t = hsi_denoised_2d_norm_t.reshape(ypix, xpix, zpix)
spectsig_plot("Four Classes of Spectral Features", hsi_denoised_3d_norm_t, labels)
# Create normalized, denoised, hyperspectral image
hsi_denoised_3d_norm = normalize(hsi_denoised_1d_transpose).reshape(ypix, xpix, zpix)
# Labeled data
n_samples = 6
n_sup = 4 # number of supervised samples
n_val = 2 # number of validation samples
hsi_X_sup = torch.stack((hsi_denoised_3d_norm_t[pts0[0][0],pts0[0][1],:],
hsi_denoised_3d_norm_t[pts0[1][0],pts0[1][1],:],
hsi_denoised_3d_norm_t[pts0[2][0],pts0[2][1],:],
hsi_denoised_3d_norm_t[pts0[3][0],pts0[3][1],:],
hsi_denoised_3d_norm_t[pts1[0][0],pts1[0][1],:],
hsi_denoised_3d_norm_t[pts1[1][0],pts1[1][1],:],
hsi_denoised_3d_norm_t[pts1[2][0],pts1[2][1],:],
hsi_denoised_3d_norm_t[pts1[3][0],pts1[3][1],:],
hsi_denoised_3d_norm_t[pts2[0][0],pts2[0][1],:],
hsi_denoised_3d_norm_t[pts2[1][0],pts2[1][1],:],
hsi_denoised_3d_norm_t[pts2[2][0],pts2[2][1],:],
hsi_denoised_3d_norm_t[pts2[3][0],pts2[3][1],:],
hsi_denoised_3d_norm_t[pts3[0][0],pts3[0][1],:],
hsi_denoised_3d_norm_t[pts3[1][0],pts3[1][1],:],
hsi_denoised_3d_norm_t[pts3[2][0],pts3[2][1],:],
hsi_denoised_3d_norm_t[pts3[3][0],pts3[3][1],:])).float()
hsi_sup_labels = torch.cat([torch.zeros(n_sup), torch.ones(n_sup),
2*torch.ones(n_sup), 3*torch.ones(n_sup)])
hsi_y_sup = pv.utils.to_onehot(hsi_sup_labels.long(), n=4) # create onehot array for supervised labels
hsi_X_val = torch.stack((hsi_denoised_3d_norm_t[pts0[4][0],pts0[4][1],:],
hsi_denoised_3d_norm_t[pts0[5][0],pts0[5][1],:],
hsi_denoised_3d_norm_t[pts1[4][0],pts1[4][1],:],
hsi_denoised_3d_norm_t[pts1[5][0],pts1[5][1],:],
hsi_denoised_3d_norm_t[pts2[4][0],pts2[4][1],:],
hsi_denoised_3d_norm_t[pts2[5][0],pts2[5][1],:],
hsi_denoised_3d_norm_t[pts3[4][0],pts3[4][1],:],
hsi_denoised_3d_norm_t[pts3[5][0],pts3[5][1],:])).float()
#hsi_X_val = torch.cat(hsi_X_val)
hsi_val_labels = torch.cat([torch.zeros(n_val), torch.ones(n_val),
2*torch.ones(n_val), 3*torch.ones(n_val)])
hsi_y_val = pv.utils.to_onehot(hsi_val_labels.long(), n=4) # create onehot array for validation labels
# Remove all labeled samples from unlabeled data by setting them to zero
hsi_X_unsup = hsi_denoised_3d_norm_t
for i in range(n_samples):
for j in range(4):
hsi_X_unsup[labels[j][i][0], labels[j][i][1], :] *= 0
hsi_X_unsup = hsi_X_unsup.reshape(xpix*ypix, zpix)
# Randomly shuffle points in image
#hsi_X_unsup = hsi_X_unsup[torch.randperm(hsi_X_unsup.size()[0])]
(hsi_loader_unsup, hsi_loader_sup, hsi_loader_val) = pv.utils.init_ssvae_dataloaders(
hsi_X_unsup, (hsi_X_sup, hsi_y_sup), (hsi_X_val, hsi_y_val), batch_size=64)
r = (len(hsi_X_sup) + len(hsi_X_val)) / (len(hsi_X_unsup))
print("Ratio of labeled data to unlabeled data: {}".format(r))
Ratio of labeled data to unlabeled data: 0.0010666666666666667
Training SSVAE Model¶
# Train the SSVAE model
hsi_ss_data_dim = (hsi_X_unsup.shape[1],)
latent_dim = 2
num_classes = 4
# Initialize model
#hsi_ssvae = pv.models.ss_reg_iVAE(
# hsi_ss_data_dim, latent_dim, num_classes, invariances=None)
hsi_ssvae = pv.models.ssiVAE(
hsi_ss_data_dim, latent_dim, num_classes, invariances=None)
# Initialize trainer
hsi_ss_trainer = pv.trainers.auxSVItrainer(hsi_ssvae, task="classification")
# Gradually increase KL weight from 1 to 2 in the first 30 epochs
hsi_ss_kl_scale = torch.linspace(1, 2, 30)
# Train model for n epochs
for e in range(50):
sc = hsi_ss_kl_scale[e] if e < len(hsi_ss_kl_scale) else hsi_ss_kl_scale[-1]
hsi_ss_trainer.step(hsi_loader_unsup, hsi_loader_sup, hsi_loader_val,
aux_loss_multiplier=50, scale_factor=sc, sampler_d="gaussian")
hsi_ss_trainer.print_statistics()
# Plot learned latent manifolds every 10 epoch
#if (e+1) % 10 == 0:
#plot_manifolds(f1_ssvae)
Epoch: 1 Training loss: 31.4974, Test accuracy: 0.2500 Epoch: 2 Training loss: 18.7537, Test accuracy: 0.2500 Epoch: 3 Training loss: 18.4341, Test accuracy: 0.2500 Epoch: 4 Training loss: 18.3721, Test accuracy: 0.2500 Epoch: 5 Training loss: 18.3608, Test accuracy: 0.2500 Epoch: 6 Training loss: 18.3342, Test accuracy: 0.2500 Epoch: 7 Training loss: 18.0592, Test accuracy: 0.2500 Epoch: 8 Training loss: 17.9969, Test accuracy: 0.2500 Epoch: 9 Training loss: 17.9891, Test accuracy: 0.2500 Epoch: 10 Training loss: 17.9900, Test accuracy: 0.2500 Epoch: 11 Training loss: 17.9861, Test accuracy: 0.2500 Epoch: 12 Training loss: 17.9681, Test accuracy: 0.2500 Epoch: 13 Training loss: 17.8745, Test accuracy: 0.2500 Epoch: 14 Training loss: 17.7099, Test accuracy: 0.2500 Epoch: 15 Training loss: 17.6864, Test accuracy: 0.2500 Epoch: 16 Training loss: 17.6844, Test accuracy: 0.2500 Epoch: 17 Training loss: 17.6963, Test accuracy: 0.2500 Epoch: 18 Training loss: 17.7018, Test accuracy: 0.2500 Epoch: 19 Training loss: 17.7058, Test accuracy: 0.2500 Epoch: 20 Training loss: 17.7134, Test accuracy: 0.2500 Epoch: 21 Training loss: 17.7141, Test accuracy: 0.2500 Epoch: 22 Training loss: 17.7331, Test accuracy: 0.2500 Epoch: 23 Training loss: 17.7405, Test accuracy: 0.2500 Epoch: 24 Training loss: 17.7440, Test accuracy: 0.2500 Epoch: 25 Training loss: 17.7537, Test accuracy: 0.2500 Epoch: 26 Training loss: 17.7608, Test accuracy: 0.2500 Epoch: 27 Training loss: 17.7585, Test accuracy: 0.2500 Epoch: 28 Training loss: 17.7826, Test accuracy: 0.2500 Epoch: 29 Training loss: 17.7857, Test accuracy: 0.2500 Epoch: 30 Training loss: 17.7846, Test accuracy: 0.2500 Epoch: 31 Training loss: 17.7967, Test accuracy: 0.2500 Epoch: 32 Training loss: 17.7903, Test accuracy: 0.2500 Epoch: 33 Training loss: 17.7868, Test accuracy: 0.2500 Epoch: 34 Training loss: 17.7914, Test accuracy: 0.2500 Epoch: 35 Training loss: 17.7877, Test accuracy: 0.2500 Epoch: 36 Training loss: 17.7859, Test accuracy: 0.2500 Epoch: 37 Training loss: 17.7848, Test accuracy: 0.2500 Epoch: 38 Training loss: 17.7857, Test accuracy: 0.2500 Epoch: 39 Training loss: 17.7896, Test accuracy: 0.2500 Epoch: 40 Training loss: 17.7852, Test accuracy: 0.2500 Epoch: 41 Training loss: 17.7861, Test accuracy: 0.2500 Epoch: 42 Training loss: 17.7836, Test accuracy: 0.2500 Epoch: 43 Training loss: 17.7865, Test accuracy: 0.2500 Epoch: 44 Training loss: 17.7867, Test accuracy: 0.2500 Epoch: 45 Training loss: 17.7840, Test accuracy: 0.2500 Epoch: 46 Training loss: 17.7890, Test accuracy: 0.2500 Epoch: 47 Training loss: 17.7753, Test accuracy: 0.2500 Epoch: 48 Training loss: 17.7842, Test accuracy: 0.2500 Epoch: 49 Training loss: 17.7806, Test accuracy: 0.2500 Epoch: 50 Training loss: 17.7772, Test accuracy: 0.2500
hsi_ssvae_z_mean, hsi_ssvae_z_sd, hsi_ssvae_z_labels = hsi_ssvae.encode(hsi_train_data)
hsi_ssvae_z_mean = normalize(hsi_ssvae_z_mean)
hsi_ssvae_z_sd = normalize(hsi_ssvae_z_sd)
hsi_ssvae_z_labels_2d = hsi_ssvae_z_labels.reshape(ypix, xpix)
hsi_ssvae_0_labels = (hsi_ssvae_z_labels_2d == 0).nonzero(as_tuple=False)
hsi_ssvae_1_labels = (hsi_ssvae_z_labels_2d == 1).nonzero(as_tuple=False)
hsi_ssvae_2_labels = (hsi_ssvae_z_labels_2d == 2).nonzero(as_tuple=False)
hsi_ssvae_3_labels = (hsi_ssvae_z_labels_2d == 3).nonzero(as_tuple=False)
Identifying Classes¶
ssvae_label_plot(hsi, hsi_ssvae_0_labels) # 180 spectra from specified class
ssvae_label_plot(hsi, hsi_ssvae_1_labels)
ssvae_label_plot(hsi, hsi_ssvae_2_labels)
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) /var/folders/2t/38t43z9j1jbd04b_hsjlptfw0000gn/T/ipykernel_92950/1746077235.py in <module> ----> 1 ssvae_label_plot(hsi, hsi_ssvae_2_labels) /var/folders/2t/38t43z9j1jbd04b_hsjlptfw0000gn/T/ipykernel_92950/2819777500.py in ssvae_label_plot(img, pts, starting_index, save) 190 191 fig_plt = fig.add_subplot(gs[row_count,0]) --> 192 fig_plt.plot(img[:, pts[i,0], pts[i,1]], 193 c=color[row_count], lw=3) 194 plt.xticks([]) IndexError: index 2 is out of bounds for dimension 0 with size 2
ssvae_label_plot(hsi, hsi_ssvae_3_labels)
Plotting Results¶
ssvae_plot("SSVAE on Artificial Data", hsi_ssvae_z_mean, hsi_ssvae_z_labels, xpix, ypix)
# X-axis
hsi_ssvae_z_mean_x = hsi_ssvae_z_mean[:,0].detach().cpu().numpy()
hsi_ssvae_z_mean_x = hsi_ssvae_z_mean_x.reshape(ypix, xpix)
print("X Min: ", (hsi_ssvae_z_mean_x==np.min(hsi_ssvae_z_mean_x)).nonzero())
print("X Max: ", (hsi_ssvae_z_mean_x==np.max(hsi_ssvae_z_mean_x)).nonzero())
hsi_ssvae_z_mean_x_sort = np.dstack(
np.unravel_index(np.argsort(hsi_ssvae_z_mean_x.ravel()), (ypix, xpix)))
hsi_ssvae_z_mean_x_sort = np.squeeze(hsi_ssvae_z_mean_x_sort)
spectra_plot("SSVAE on Artificial Data X-Axis Spectra", hsi, hsi_ssvae_z_mean_x_sort, xpix, ypix)
X Min: (array([25, 25]), array([53, 54])) X Max: (array([ 0, 0, 0, ..., 149, 149, 149]), array([ 0, 1, 2, ..., 147, 148, 149]))
# Y-axis
hsi_ssvae_z_mean_y = hsi_ssvae_z_mean[:,1].detach().cpu().numpy()
hsi_ssvae_z_mean_y = hsi_ssvae_z_mean_y.reshape(ypix, xpix)
print("Y Min: ", (hsi_ssvae_z_mean_y==np.min(hsi_ssvae_z_mean_y)).nonzero())
print("Y Max: ", (hsi_ssvae_z_mean_y==np.max(hsi_ssvae_z_mean_y)).nonzero())
hsi_ssvae_z_mean_y_sort = np.dstack(
np.unravel_index(np.argsort(hsi_ssvae_z_mean_y.ravel()), (ypix, xpix)))
hsi_ssvae_z_mean_y_sort = np.squeeze(hsi_ssvae_z_mean_y_sort)
spectra_plot("SSVAE on Artificial Data Y-Axis Spectra", hsi, hsi_ssvae_z_mean_y_sort, xpix, ypix)
Y Min: (array([35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37,
38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41,
41, 41, 41, 41, 41, 42, 42, 42, 42, 43, 43, 44, 44, 56, 56, 58, 58,
59, 59, 61, 61, 61, 64, 64, 64, 69, 69, 73, 74, 74, 74, 77, 77, 79,
79, 80, 80, 80, 80, 80, 80, 80, 80, 80, 81, 81, 81]), array([75, 76, 77, 78, 79, 71, 72, 73, 74, 85, 86, 87, 67, 68, 69, 70, 77,
62, 63, 64, 65, 66, 85, 87, 58, 59, 60, 71, 55, 56, 57, 93, 94, 49,
50, 51, 52, 95, 96, 58, 85, 98, 99, 50, 51, 84, 85, 48, 49, 54, 55,
55, 56, 59, 60, 61, 64, 65, 66, 69, 70, 77, 73, 74, 75, 75, 76, 92,
93, 80, 81, 82, 83, 86, 87, 88, 89, 90, 90, 91, 92]))
Y Max: (array([ 40, 40, 111, 111]), array([14, 15, 73, 74]))
custom_spect_labels = [hsi_ssvae_0_labels[0], hsi_ssvae_1_labels[0],
hsi_ssvae_2_labels[1], hsi_ssvae_3_labels[0]]
# Plot one spectrum from each class in the latent representation of the data
custom_spect_label_plot("Classes from SSVAE on Artificial Data", hsi, custom_spect_labels)